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NEW RAPID ALGORITHMS FOR FACTOR ANALYSIS BY UNWEIGHTED LEAST SQUARES, 



GENERALIZED LEAST SQUARES AND MAXIMUM LIKELIHOOD 
Karl G. Jtireskog and Marielle van Thillo 

1 . Introduction 

We shall describe a new basic algorithm that may be used to do factor 
analysis by any of the three methods 

(i) unweighted least squares (ULS) 

(ii) generalized least squares (GI£>) 

(iii) maximum likelihood (ML). 

The ULS method produces solutions that are equivalent to those obtained 
by the iterated principal factor method and the minres method (see Harman, 

1967, Chapters 8 and 9 )* The generalized least squares method is described 
by J&reskog and Goldberger (1971)* In the ML case, the new algorithm is 
simpler and faster than J&reskog’ s (1967a, b) method UMLFA. It is similar 
to Clarke’s (1970) algorithm but Heywood cases are handled in a simpler and 
more efficient way. Although the new algorithm handles ULS and GLS as 
well as ML, the computer program is shorter than UMLFA. 

The GLS and ML methods are scale free. When multivariate normality is 
assumed, both GIE and ML yield estimates that are asymptotically efficient. 
Both GIE and ML require a positive definite variance-covariance matrix S 
or correlation matrix R ; UIE will work even on a matrix that is non-Gramian. 

The model is the usual factor analysis model, which requires the popu- 
lation variance -covariance matrix or correlation matrix £ of the observed 
variables to be of the form 
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E = AA' + y , (l) 

2 

where A is a p x k matrix of factor loadings and is a p x p 

diagonal matrix of unique variances. The factors are assumed to be 
orthogonal. 

The model (l) is fitted to the observed variance -covariance matrix 
S or to the corresponding correlation matrix R , by the minimization 
of a fitting function F(A,t|r) , which is different for each of the three 
methods. The minimization of. F(A, \|r) is done in two steps. First the 
conditional minimum of F for given \|r is found. This, gives a func- 
tion f(\|r) which is then minimized numerically using the Newton-Raphson 
procedure. Function values and derivatives of f of first and second 
order are given in terms of the characteristic roots and vectors of a 

certain matrix A . In the GLS and ML cases a transformation from t. to 

1 

0. is made to obtain stable derivatives at ty. = 0 . The basic formulas 

l l 

for ULS, GIB and ML are given in sections 2 , 3 and 4 respectively. 
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2. Formulas for ULS 
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Fitting function: F(A,\|/) = (l/2) tr[ (S - Z) ] , 

Z = AA* + f 2 , 

A 1 A is assumed to be diagonal 

r 

Matrix whose roots and vectors are computed: A « S - \|r 



Characteristic roots: y > y > 



> 






.1/2 



Corresponding orthonormal vectors: * * ‘ >l p 

Conditional solution for A for given \|r : A = fl I 1 ^ ** } 

r 1 = diag( 7;L ,7 2 , . 
Of = CojjOJg . . .ox ] 

Function minimized by the Newton-Raphson method: 



( 2 ) 

(3) 

oo 



(5) 



•, 7 V ) , 



f(t) “ (1/2) 2 7 



m=k+l 



m 



( 6 ) 



First order derivatives: bf/ht. = -2^. Z 7 m . 

1 1 , , . 'm 1 

m=k+l 



2 

f • 

im 



Second order derivatives: 



( 7 ) 



d'^f/dt.d*. = btf.f. Z ©. oj. Z 
' 1 .1 in , . , im im 



k 7 m + 7 



m ' n 



CD. CD. + 8 . 



1 , , , ira jra .7-7 in in ii 

0 m=k+l 0 n=l 'm 'n 0 m=k+l 



(#? - 7 / 2 )cd? ] ( 8 ) 

v 1 7 m' 7 im J x 7 



P Jr 2 

Approximate second order derivatives: 8 f/S^.S^. = Z co. co. ) (9) 

1 d 1 d m=k+l d 
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3 • Formulas for GLS 



Fitting function: F(A, \|r) = (l/2) tr[(S ^Z - i) 2 ] , 

Z s= AA' + \|r , 

_2 

A'\|r A is assumed to be diagonal 
Matrix whose roots and vectors are computed: A = \|/S \|r 



Characteristic roots: < y n < ... <7 

' 1 — ' 2 — — / 



P 



Corresponding orthonormal vectors: . 

Conditional solution for A for given \|r : 



,cd 



P 



a = ^( r " 1 - i ) 1 / 2 , 

r x = diag(7 x ,7 2 , *“>7 k ) 



= [oo^g. . .0^] 

nr „ 

Transformation: f. = + ye 1 ’ = ■*-°® 1 

Function minimized by the Newton-Raphson method: 



f( 0 ) = (1/2) Z ( 7 m - l ) 2 



m=k+l 



First order derivatives: bf/bo^ = 



m=k+l 



(7 



m 



7 )cd? 
'm' 1 m 



Second order derivatives: 



_ -n k 7 + 7 - 2 

.2 „ /N _ N ~ £ r £ 'm 'n 



b f/bd.bd. = B. .bf/bQ. + Z 7 CD. CD. [ Z 7 
' in id 1 . 'ra ira Dm , n 

u 0 m=k+l 0 n=l 



CD. CD . + 



7-7 m D n 

'm ' n 0 



2 

Approximate second order derivatives: 8 f /be be. = 

J- t J 



P P 

( Z CD. CD. ) 
m=k+l im 



lERJC 



(10) 

(11) 

(12) 

(13) 

(14) 

(15) 

(16) 

1 D 

(18) 
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4 . Formulas for ML 



Fitting function: F(A,\|r) = tr(E _1 S) - log | S _J_ S | - p , 






(19) 



E = AA' + t 



_2 

A’ \|r A is assumed to be diagonal 



Matrix of which roots and vectors are computed: A = \JrS \|/ 



Characteristic roots: < ... <y 

' 1 — ' 2 — — X 



Corresponding ortho normal vectors: 
Conditional solution for A for given \|r 



Transformation: * , 0^ = log ^ 

Function minimized by the Newton -Raphson method: 



.CO 



a = t^Cr" 1 - i) 1 / 2 



r 1 = diag( 7;L ,7 2 , . . .,/ k ) 



ft, = [cOjCOg. . .01 ] 



f(0) = s (log 7rn + 1/7 - 1) 

. , , 'm ' 'm 

m=k+l 

P 2 

First order derivatives: Sf/S0. = E (l - 1/7 )co. 

' 1 , , v ''m' 1m 

m=k+l 



( 20 ) 

( 21 ) 



Second order derivatives: 



o . p 

d t/be.be. = -B_. .Sf/Se. + E co. co. [ E 
' 1 n xj 1 , . , im nm 



k ^m 4 ^n 2 c t 

co. 00 . + 5. J 

m=k+l a® 7 m - 7 n in jn xj 



( 22 ) 



Approximate second order derivatives: S 2 f/Se.Se ± ( E co co ) 2 ( 2 7 >) 

1 J m=k+l im J‘ m 
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V 


5* Basic Minimization Algorithm 

Let 0 denote a column vector with elements 0^, 0^, • • • , 0p (GLS and 

ML) or f^ } f } • • • > & (ULS), and let h and H denote the column vector 

2 

and matrix of corresponding derivatives dg/d© and 5 g /<5q<5q' , respectively. 

( s ) ( s ) * 

Let 0' ' denote the value of 0 in the s iteration and let h v and 

(s) 

H' ' be the corresponding vector and matrix of first- and second-order , 

derivatives. The iteration procedure may then be written 

H (s) s( s > = h (s) , (24) 

0 (s+l) = 0 (s) _ 8 (s) ^ (25) 

where &t s ) is a column vector of corrections determined by (24). The 
Newton-Raphson procedure is therefore easy to apply, the main computations 
in each iteration being the computation of the roots and vectors of A 
and the solution of the symmetric system (24). It has been found that the 
Newton-Raphson procedure is very efficient, generally requiring only a few 
iterations for convergence. The convergence criterion is that the largest 
absolute correction be less than a prescribed small number e . The 
minimizing 0 may be determined very accurately, if desired, by choosing 
6 very small. 


r 

t. 

> 

i. 

i 


In detail, the numerical method is as follows : the starting point 

0^^ is chosen as (see e.g., jBreskog, 1963, eqs. 6.20 and 7*10 or 


L 

<■ 


JSreskog, 1967* eq. 26), 

a* 1 ) = io g [(i - k/apj/s 11 ] , = +J[(i - k/apj/s 11 ] , 


(: 

V- 

t-i 

f 

?■ 

p 


ft 1 ) = ,6s.. if S is not (26) 

i li 

positive definite, 


1 
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positive definite, 
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where s^ and s^ are the i^ 1 diagonal element of S and S ^ respec- 
tively. The exact matrix H of second order derivatives given by (8), (17) 
or (22) may not be positive definite in the beginning. Therefore, the 
approximation E given by (9)., (l8) or (23) is used in the first iteration 
and for as long as the maximum absolute correction is greater than a given 
constant e (see sec. 7)* After that, H is used if it is positive 

ill 

definite. It has been found empirically that E gives good reductions in 
function values in the early iterations but is comparatively ineffective 
near the minimum, whereas H near the minimum is very effective. 

In each iteration we compute the characteristic roots and vectors of 
A by the Householder transformation to tridiagonal form, the OR method for 
the roots of the tridiagonal matrix and inverse iteration for the vectors. 
This is probably the most efficient method available (see Wilkinson, 1965)- 
The system of equations (24) is solved by the square root factorization 
H = TT' , where T is lower triangular. This shows at an early stage 
whether H is positive definite or not. 

In Heywood cases, when one or more of the 0^ -» _0 ° > i.e., 1*4 -» 0 , 

a slight modification of the Newton-Raphson procedure is necessary to achieve 
fast convergence. This is due to the fact that the search for the minimum 
is then along a "valley" and not in a quadratic region. For ML and GLS, 

p 

when e. -»-<», c)g/c)0. -» 0 and 8 g/80.80. 0 , j = 1, 2, ...,p , 

so that when 0. is small the i^ element of h and the i^* 1 row and 
1 

column of H and E are also small. This tends to produce a "bad" 
correction vector S and the function may increase instead of decrease. 

A simple and effective way to deal with this problem is to delete the i^* 1 



O 
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equation in the system (24) and compute the corrections for all the other 
0's from the reduced system. One then computes the correction for 0_^ 
as 

8 i = (Sg/^0 i )/(S 2 g/Se?) • (27) 

This procedure will decrease 0^ slowly in the beginning but faster 

the more evident it is that 0. is a Heywood variable. When 0. has 

x 1 

become less than log(e) it is not necessary to change 0_^ any more unless 
5f/50^ is negative. For UI£, an analogous procedure is used. When \|r. 
becomes less than , \|n is not changed unless Bf/d\|/^ is negative. 

Thus, the procedure corrects itself quickly if a variable is incorrectly 
taken as a Heywood variable. 
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6. The Program 



In this section we describe briefly what the program does- Details 
about the input are given in section 7* For those users who feel too re- 
stricted in their choice of an input matrix, as provided by the program, the 
kernel of the program is available as a subroutine. The input and output 
parameters for that subroutine will be described in section 8. 

The input data may be raw data from which the matrix to be analyzed 
is computed, or it may be a dispersion matrix, or it may be a correlation 
matrix or a correlation matrix followed by a vector of standard deviations. 
From these input matrices, variables may be selected to be included in the 
analysis, so that the matrices to be analyzed could be of smaller order than 
the input matrices. Variables may also be interchanged with one another. 

The matrices to be analyzed may be dispersion matrices or correlation 
matrices. The user has the option to read in a starting point for \|r or 
have the program define a starting point (see sec. 5 )* This can be useful 
if convergence is slow and the user runs out of computer time- From the 
intermediate results the last \|/ can be read in as a new starting point 
and minimization can continue. 

For the given matrix S to be analyzed of order p by p and a 
given lower bound k^ and a 

factors, the program performs a sequence of factor analyses by the ML, 

ULS or GLS method of estimation chosen by the user and outlined in the 
previous sections. One such analysis is done for each number of factors 

k = k^, k^+1, • • • , ky 



given upper bound k^ 



iwx one iiUIUUGX 
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The output will consist of the title with parameter listing and the 
matrix to be analyzed. Then for each number of factors k the unrotated 

factor loadings, the unique variances and the varimax-rotated factor loadings 

. 2 
are printed. For ML and GLS this is followed by X and the corresponding 

K. 

degrees of freedom d^ , the probability level, i.e., the probability of 

2 

obtaining a larger value of X than that actually obtained given that the 
model and the assumptions hold, and Tucker and Lewis' (1970) reliability 
coefficient , defined as follows 

C Q = N - 1 (l/ 6 )( 2 p + 5 ) 

2 P . , 

X Q = C 0 C E log s i;L - log | S | ] 

i=l 

d o = I p ( p ' ^ 

“o ■ x o/ a o 
C k = C 0 - (2/3)* 

X 2 = C v f • 
k k min 

d k = | [ ( p " k) 2 - (p + k)] 

D M 0 ^ 

k = M 0 - ! 

Finally the latent roots and their first differences at the minimum and the 
matrix of residual correlations are printed. The user also has an option to 







) 
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print intermediate results consisting of the value of the function and the 
vector \|r at each iteration. Examples of input and output can be found 
in Appendix B. 

The following limitations are imposed on the program: 
max. no. of variables after selection = 50 
max. no. of variables before selection = 75 
max. no. of factors = 50 

storage requirements on the IBM 560/65 = 120K ( K = 1024 bytes ) 

The program can easily be modified to allow for a larger number of variables 
and factors. Instructions on how to change the maximum number of variables 
and factors allowed by the program can be found in Appendix C. The pro- 
gram is written in FORTRAN IV -G and has been tested out on the IBM 560/65 
at Educational Testing Service. Double precision is used in floating point 
arithmetic throughout the program. With minor changes the program should 
run on any computer with a FORTRAN IV compiler. In computers with a single 
word length of 5 6 bits or more, single precision is probably sufficient. 

Although the program has been working satisfactorily for all data 
analyzed so far, no claim is made that it is free of error and no warranty 
is given as to the accuracy and functioning of the program. 




bl 
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7 . Input Data 

For each set of data to be analyzed, the input consists of the following: 

a. Title card 

b. ' Parameter card 

c . Data matrix 

d. Selection cards (optional) 

e. Starting point (optional) 

f . New data or a STOP card 

The function and setup of each of the above quantities are described in 
general terms below. Illustrative examples are given in Appendix B. 

a. Title Card 

Whatever appears on this card will appear on the first page of the 
printed output. All 80 columns of the card are available to the user. 

b. Parameter Card 

All quantities on this card, except for the logical indicators, must be 
punched as integers right adjusted within the field, 
cols. 1-5 
cols. 6-10 
cols. 11-15 
cols. 16-20 
cols. 21-25 

col. 51 



number of observations N 

order of data matrix ( p^ ), before selection of variables 

lower bound for the number of factors k^ 

upper bound for the number of factors k^ 

maximum number of iterations allowed for each number 

of factors k 

logical variable which determines whether selection 
of variables from the data matrix is desired 
col. 31 = T, if selection of variables is wanted 
col. 31 = F, if no selection of variables is wanted 
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col. 



col. 



col. 







32 logical variable which determines whether a dis- 

persion matrix or a correlation matrix is to be 
analyzed 

col. 32 = T, if a dispersion matrix is to be 
analyzed 

col. 32 = F, if a correlation matrix is to be 
analyzed 

41 integer indicator which determines whether raw data, 

a dispersion matrix, a correlation matrix or a 
correlation matrix with standard deviations are 
read in to determine the matrix to be analyzed 

col. hi = 1, read in raw data 

col. 4l = 2, read in a dispersion matrix 

col. 4l = 3, read in a correlation matrix 

(followed by a vector of standard 
deviations if col. 32 is T) 

42 integer indicator which determines which method 
of estimation is to be used 
col. 42 = 1 for ULS 
col. 42 = 2 for GIS 
col. 42 = 3 for ML 

integer indicator which determines whether inter- 
mediate results are to be printed 
col. 43 = 0, if no_ intermediate results are to be 
printed 

col. 43 = 1, if intermediate results are to be 
printed (see sec. 6) 



col. 43 
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col. 44 



cols. 46-55 



cols. 56-65 



cols. 66-70 



integer indicator which determines whether a 
starting point is defined by the program or is to 
be supplied by the user (see secs. 5 and 6 ) 
col. 44 = 0 , if a starting point. is defined by the 
program 

col. 44 = 1 , if a starting point is read in as data 

convergence criterion e (see sec. 5)* For 

reasonable results use e < .005 • 

e_ : if all elements of the correction vector are 
E 

less than e_ the exact second order derivatives 

E 

are computed in the minimization algorithm, other- 
wise the approximate second order derivatives are 
used (see sec. 5)* From our experience e_ = .1 
seems reasonable. 

logical tape (disk) number of scratch tape (disk) 
used for intermediate storage 



c. Data Matrix 

The data matrix is preceded by a format card, containing at most 80 
columns, beginning with a left parenthesis and ending with a right paren- 
thesis. The format must specify floating point numbers consistent with the 
way in which the elements of the matrix are punched. Users who are un- 
familiar with FORTRAN are referred to a FORTRAN Manual where format rules 



are given. 
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The input matrix can be any one of the following: 

If col. In a 1 on the parameter card an N x p matrix of raw data is 
read in, one row at a time, starting a new card for each row. The matrix 
is preceded by a format card as described above. 

If col. 4l = 2 the lower triangular part of a dispersion matrix, 
including the diagonal, is read in. The matrix should be punched row -wise 
as one long vector, i.e., there is no need to go to a new card if a new row 
starts. Again the matrix should be preceded by a format card. 

If col. 4l = 3 and col. $2 = F the lower triangular part of a cor- 
relation matrix, including the diagonal, is read in. The matrix should be 
punched row-wise as one long vector, and should be preceded by a format card. 

If col. 4l = 3 and col- $2 = T the lower triangular part of a cor- 
relation matrix, including the diagonal, is read in. The matrix should be 
punched row-wise as one long vector, and should be preceded by a format card. 
This matrix is then followed by a format card and a row-vector of standard 
deviations . 

d. Selection Cards (optional) 

Omit if column 31 of the parameter card is F. Otherwise the first card 
will have an integer value p punched in columns 1-5 > right adjusted 
within the field. This integer will specify the new order of the data 
matrix after selection of variables ( P < p^ )• 

The next card will contain integers, right adjusted in five column 
fields (i.e., sixteen such values will fit on one card), specifying which 
columns (rows) are to be included . For example, if ]n = 6 , p = 3 and 
the first, second and fifth columns (rows) are to be excluded, this card 




1 6 



would have a 3 punched in column 5j a 4 punched in column 10 and a 6 punched 
in column 15 • 

Note that if p = p Q there will be no reduction in the size of the 
data matrix but columns (rows) can be interchanged. 

e . Starting Point (optional) 

Omit if column 44 on the parameter card is zero. Otherwise read in a 
starting \j/ vector punched according to a format of 5 D 15 «T (see sec. 6) 
for each number of factors k . 

f. Stacked Data 

In the preceding paragraphs we have described how each set of data 
should be set up. Any number of such sets of data may be stacked together 
and analyzed in one run. After the last set of data in the stack, there 
must be a card with the word STOP punched in columns 1-4. 
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8. Subroutine NWTRAP 

As an alternative to the program the user can write his own main 
program in which he calls the minimization package with a CALL to NWTRAP. 
The following subroutines are called by NWTRAP and thus are part of the 
package: FCTGR, INCPSI, S0LVE, ISMSL, HFWLIN, TRIDI, EIGVEC and QRB. 

The minimization package should be called with the following sequence 
of F0RTRAN statements: 

D0 10 K = KL,KU 

CALL NWTRAP (P, K, II, 12, 13, S, EPS, EPSE, MAXIT, A, E, X, Y, FO, DET ) 

• 

10 CONTINUE 

The D0 loop runs from KL, tho lower bound on the number of factors to KU, 
the upper bound on the number of factors. 

Next follows a description of the parameters of NWTRAP: 

Input Parameters 

P order of data matrix S 

K number of factors 

II determines which method of estimation is to be used 



11 = 


1 


for 


ULS 


11 = 


2 


for 


GLS 


11 = 


3 


for 


ML 



12 determines whether intermediate results are to be 

printed (see sec. 6) 

12 = 0 , if no intermediate results are to be printed 
12 = 1 , if intermediate results are to be printed 
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13 


determines whether the starting point is defined 
by the program or read in by the user as data 
(see sec. 6) 

13 = 0 , if the starting point is defined by the 
program (see sec. 5) 

13 - 1 ) a starting vector t is read in as data 
with a format of 5D15*7 for each number 
of factors K 


S 


data matrix, stored row-wise as a vector. Should 
be singly dimensioned in the calling program by at 
least (P(P + l))/2 . 


EPS 


convergence criterion e (see sec. 5). For reason 
able results use e < .005 • 


EPSE 


if all elements of the correction vector are less 
than EPSE, the exact second order derivatives are 
used in the minimization algorithm, otherwise the 
approximate second order derivatives are used. 
From our experience EPSE = .1 seems reasonable. 


MAXIT 


maximum number of iterations allowed for each 
number of factors K . Program exits if this 
number is exceeded. 


Output Parameters 





A 



matrix of unrotated factor loadings, stored row- 
wise as a vector. Should be singly subscripted 
in the calling program by at least P x K . 
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E dummy vector. Should be dimensioned in the calling 

program by at least (P x (P + l))/2 . 

X vector of unique variances. Should be dimensioned 

in the calling program by at least P . 

Y vector of latent roots of A at the minimum. 

Should be dimensioned in the calling program by at 
least P • 

FO the value of the function at the minimum 

DET the determinant of the data matrix S 
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Appendix B 



An Example of Input Data 

We shall illustrate how input data should be set up by means of four 
sets of data. In all four cases Harman's correlation matrix of twenty-four 
psychological tests for 145 children is the input matrix. In the first set 
of data all twenty-four variables are to be analyzed with 4 and 5 common 
factors using the ML method of estimation. Intermediate output is requested. 

The second set of data is as the first except that the ULS method of 
estimation is used and no intermediate output is to be printed. 

The third set of data selects the first 15 variables from the input 
matrix to be analyzed with 4 common factors and using the ML method of 
estimation. No intermediate output is to be printed. 

The last set of data analyzes the matrix of the third data set with 
the Heywood variable (the 11th variable) removed. Thus 12 variables are 
selected from the first thirteen and are analyzed with 5 and 4 common 
factors using the ML method of estimation. The 12 variables selected are 
also rearranged so that the variables appear in the following order: 

15, 10, 12, 5, 4, 9, l, 8, 7, 2, 6, 5 

i.e., the 15th variable is now the first, the 10th variable is now the 
second, etc. 

The next two pages show card by card how the data should be punched. 

One line corresponds to one card. For all sets of data MAXIT has been 
set to 50 and the scratch unit is 4. 

The results obtained from these data follow on subsequent pages. 
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Appendix C 

If the user wishes to change the number of variables, p , and/or the 
number of variables before selection, , the MAIN program and subroutines 

REX, SELECT, NWTRAP, INCPSI and FCTGR need to be modified. 

In the MAIN program the DIMENSI0N card should read as follows: 

DIMENS I0N FMT (lO),S(n),A(n), HEAD (20),YY(p),E(m),Y(p) 
where n = (p Q (p 0 + l))/2 and m = (p(p + l))/2 . 

In subroutine REX the DIMENSI0N card should be: 

DIMENS I0N S(l),E(l),Y(p 0 ),X(p 0 ),FMT(lO) . 

In subroutine SELECT the DIMENSI0N card should be: 

DIMENS I0N S(l),E(l),MM(p 0 ) • 

In subroutines NWTRAP, INCPSI and FCTGR the C0MM0N block KERN should 

read: 

C0MM0N/KERN/G(p), V(p), VB(p), D2(p), Sl(p), S2(p), S3 (p),EPSU, 

BND, IM(p ) , M0R, KP, M0RE, MAXTRY, P2, KP1 . 



Caution : The following relationship between p and p^ must hold 

P 2 5 (P 0 ( p o + 1 ))/ 2 
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